function est = fitfsmod(D,modspecs)

% -- Add intercept to data if not included--;
if ~any(strcmp(D.Properties.VariableNames,'intercept'))
    D.intercept(:) = 1;
end

% -- Check to see if all variables in model are in dataset and non missing --;
modvars = [modspecs.outcome modspecs.controls modspecs.fsfevars];
if ~all(ismember(modvars,D.Properties.VariableNames))
    error('Variables in modspecs not present in data')
elseif any(isnan(reshape(D{:,modvars},[],1)))
    error('Variables in model contain missing values in data')
elseif ~all(ismember(modspecs.fsfevars,modspecs.controls))
    error('Teacher fixed effect variables also need to be listed in controls')
end

% -- If Random Effects Model --;
if isempty(modspecs.fsfevars)
    fsmdl = fitlm(D,strcat(modspecs.outcome,'~',strjoin(modspecs.controls,'+')),'Intercept',false);
    controls = renamevars(fsmdl.Coefficients,{'Estimate','SE'},{'Coefficient','StandardError'});
    est = struct('outcome',modspecs.outcome...
        ,'n',fsmdl.NumObservations...
        ,'numparms',fsmdl.NumCoefficients...
        ,'controls',controls(modspecs.controls,{'Coefficient','StandardError'})...
        ,'residvar',fsmdl.MSE);
    return
end

% -- Outcomes --;
n = size(D,1);

Y = D.(modspecs.outcome);

% -- Controls --;
control_varnames = modspecs.controls(~ismember(modspecs.controls,modspecs.fsfevars));
X = D{:,control_varnames};

% -- Teacher Effect Variables --;
fevars = modspecs.fsfevars;
num_fe = numel(fevars);
Z = +D{:,fevars}; 
felvl = round((Z'*Z) \ (Z'*ones([size(Z,1) 1])),10)'==1;
if ~all(sum(Z(:,felvl),2)==1)
    error('Level variables cannot be identified')
end

[uni_teacherid,~,tch_id] = unique(D.teacherid);
numTCH = numel(uni_teacherid);

feIND = arrayfun(@(j) fullcolrank(Z(tch_id==j,:)),1:numTCH,'unif',0);
feIND = cat(1,feIND{:});

Z = arrayfun(@(c) sparse((1:n),tch_id,Z(:,c),n,numTCH),1:num_fe,'unif',0);
Z = cat(2,Z{:});

Z = Z(:,feIND(:));

% -- Estimation --;

ZX = Z'*X;
ZZinv = (Z'*Z) \ speye(size(Z,2));

MZ = speye(numel(Y)) - Z*ZZinv*Z';
b = (X'*MZ*X) \ (X'*MZ*Y);

fe = ZZinv*(Z'*(Y - X*b));

% -- standard Errors --;
resid = (Y - X*b - Z*fe);

num_parms = numel(b) + numel(fe);
residvar = sum(resid.^2)/(n - num_parms);

KM = (X'*X - ZX'*(ZZinv*ZX)) \ eye(numel(b));

b_var = residvar*diag(KM);

% -- Output --;
est = struct('outcome',modspecs.outcome,'n',n,'numparms',num_parms...
    ,'controls',[]...
    ,'residvar',residvar);
est.controls = array2table(nan([numel(modspecs.controls) 2])...
    ,'VariableNames',{'Coefficient','StandardError'}...
    ,'RowNames',modspecs.controls);
est.controls{control_varnames,'Coefficient'} = b;
est.controls{fevars,'Coefficient'} = 0;
est.controls{control_varnames,'StandardError'} = sqrt(b_var);
